R Notebook

This is an R Markdown Notebook. When you execute code within the notebook, the results appear beneath the code.

Try executing this chunk by clicking the Run button within the chunk or by placing your cursor inside it and pressing Ctrl+Shift+Enter.

# example code from textbook
library(astsa)
par(mfrow=2:1)
tsplot(jj, ylab="QEPS", type="o", col=4, main="Johnson & Johnson

Quarterly Earnings")

tsplot(log(jj), ylab="log(QEPS)", type="o", col=4)

library(xts)
djia_return = diff(log(djia$Close))[-1]
par(mfrow=2:1)
plot(djia$Close, col=4)
plot(djia_return, col=4)

tsplot(diff(log(gdp)), type="o", col=4, ylab="GDP Growth") # diff-log
points(diff(gdp)/lag(gdp,-1), pch=3, col=2)

par(mfrow = c(2,1))

tsplot(soi, ylab="", xlab="", main="Southern Oscillation Index", col=4)
text(1970, .91, "COOL", col="cyan4")

text(1970,-.91, "WARM", col="darkmagenta")
tsplot(rec, ylab="", main="Recruitment", col=4)

culer = c(rgb(.85,.30,.12,.6), rgb(.12,.67,.86,.6))
tsplot(Hare, col = culer[1], lwd=2, type="o", pch=0,
ylab=expression(Number~~~(""%*% 1000)))
lines(Lynx, col=culer[2], lwd=2, type="o", pch=2)

legend("topright", col=culer, lty=1, lwd=2, pch=c(0,2),

legend=c("Hare", "Lynx"), bty="n")

par(mfrow=c(3,1))
culer = c(rgb(.12,.67,.85,.7), rgb(.67,.12,.85,.7))
u = rep(c(rep(.6,16), rep(-.6,16)), 4) # stimulus signal
tsplot(fmri1[,4], ylab="BOLD", xlab="", main="Cortex", col=culer[1],

ylim=c(-.6,.6), lwd=2)
lines(fmri1[,5], col=culer[2], lwd=2)

lines(u, type="s")
tsplot(fmri1[,6], ylab="BOLD", xlab="", main="Thalamus", col=culer[1],

ylim=c(-.6,.6), lwd=2)

lines(fmri1[,7], col=culer[2], lwd=2)
lines(u, type="s")

tsplot(fmri1[,8], ylab="BOLD", xlab="", main="Cerebellum",
col=culer[1], ylim=c(-.6,.6), lwd=2)
lines(fmri1[,9], col=culer[2], lwd=2)

lines(u, type="s")
mtext("Time (1 pt = 2 sec)", side=1, line=1.75)

par(mfrow=2:1)
w = rnorm(500) # 500 N(0,1) variates
v = filter(w, sides=2, filter=rep(1/3,3)) # moving average
tsplot(w, col=4, main="white noise")
tsplot(v, ylim=c(-3,3), col=4, main="moving average")

set.seed(90210)
w = rnorm(250 + 50) # 50 extra to avoid startup problems
x = filter(w, filter=c(1.5,-.75), method="recursive")[-(1:50)]
tsplot(x, main="autoregression", col=4)

set.seed(314159265) # so you can reproduce the results
w = rnorm(200); x = cumsum(w) # random walk
wd = w +.3; xd = cumsum(wd) # random walk with drift
tsplot(xd, ylim=c(-2,80), main="random walk", ylab="", col=4)
abline(a=0, b=.3, lty=2, col=4) # plot drift

lines(x, col="darkred")
abline(h=0, col="darkred", lty=2)

t = 1:500
cs = 2*cos(2*pi*(t+15)/50) # signal
w = rnorm(500) # noise
par(mfrow=c(3,1))
tsplot(cs, col=4, main=expression(2*cos(2*pi*(t+15)/50)))
tsplot(cs+w, col=4, main=expression(2*cos(2*pi*(t+15)/50+N(0,1))))

tsplot(cs+5*w,col=4, main=expression(2*cos(2*pi*(t+15)/50)+N(0,5^2)))

ACF = c(0,0,0,1,2,3,2,1,0,0,0)/3
LAG = -5:5
tsplot(LAG, ACF, type="h", lwd=3, xlab="LAG")

abline(h=0)

points(LAG[-(4:8)], ACF[-(4:8)], pch=20)
axis(1, at=seq(-5, 5, by=2))

x = rnorm(100)
y = lag(x,-5) + rnorm(100)
ccf(y, x, ylab="CCovF", type="covariance", panel.first=Grid())

(r = acf1(soi, 6, plot=FALSE)) # sample acf values

[1]  0.60410089  0.37379533  0.21412447  0.05013659 -0.10703704 -0.18698742

par(mfrow=c(1,2), mar=c(2.5,2.5,0,0)+.5, mgp=c(1.6,.6,0))
plot(lag(soi,-1), soi, col="dodgerblue3", panel.first=Grid())
legend("topleft", legend=r[1], bg="white", adj=.45, cex = 0.85)

plot(lag(soi,-6), soi, col="dodgerblue3", panel.first=Grid())
legend("topleft", legend=r[6], bg="white", adj=.25, cex = 0.8)

set.seed(101011)
x1 = sample(c(-2,2), 11, replace=TRUE) # simulated coin tosses
x2 = sample(c(-2,2), 101, replace=TRUE)
y1 = 5 + filter(x1, sides=1, filter=c(1,-.5))[-1]
y2 = 5 + filter(x2, sides=1, filter=c(1,-.5))[-1]
tsplot(y1, type="s", col=4, xaxt="n", yaxt="n") # y2 not shown
axis(1, 1:10); axis(2, seq(2,8,2), las=1)

points(y1, pch=21, cex=1.1, bg=6)

acf(y1, lag.max=4, plot=FALSE)

Autocorrelations of series ‘y1’, by lag

     0      1      2      3      4 
 1.000 -0.223 -0.623  0.219  0.416 

acf(y2, lag.max=4, plot=FALSE)

Autocorrelations of series ‘y2’, by lag

     0      1      2      3      4 
 1.000 -0.435  0.043  0.061 -0.075 

par(mfrow=c(3,1))
acf1(soi, 48, main="Southern Oscillation Index")

 [1]  0.60  0.37  0.21  0.05 -0.11 -0.19 -0.18 -0.10  0.05  0.22  0.36  0.41  0.31  0.10 -0.06 -0.17 -0.29 -0.37 -0.32 -0.19 -0.04
[22]  0.15  0.31  0.35  0.25  0.10 -0.03 -0.16 -0.28 -0.37 -0.32 -0.16 -0.02  0.17  0.33  0.39  0.30  0.16  0.00 -0.13 -0.24 -0.27
[43] -0.25 -0.13  0.06  0.21  0.38  0.40

acf1(rec, 48, main="Recruitment")

 [1]  0.92  0.78  0.63  0.48  0.36  0.26  0.18  0.13  0.09  0.07  0.06  0.02 -0.04 -0.12 -0.19 -0.24 -0.27 -0.27 -0.24 -0.19 -0.11
[22] -0.03  0.03  0.06  0.06  0.02 -0.02 -0.06 -0.09 -0.12 -0.13 -0.11 -0.05  0.02  0.08  0.12  0.10  0.06  0.01 -0.02 -0.03 -0.03
[43] -0.02  0.01  0.06  0.12  0.17  0.20

ccf2(soi, rec, 48, main="SOI & Recruitment")

set.seed(1492)
num = 120
t = 1:num
X = ts( 2*cos(2*pi*t/12) + rnorm(num), freq=12 )
Y = ts( 2*cos(2*pi*(t+5)/12) + rnorm(num), freq=12 )
Yw = resid(lm(Y~ cos(2*pi*t/12) + sin(2*pi*t/12), na.action=NULL))
par(mfrow=c(3,2))
tsplot(X, col=4); tsplot(Y, col=4)
acf1(X, 48); acf1(Y, 48)

 [1]  0.51  0.32 -0.03 -0.31 -0.58 -0.59 -0.52 -0.25  0.02  0.36  0.46  0.64  0.51  0.26 -0.09 -0.30 -0.53 -0.58 -0.45 -0.21  0.03
[22]  0.36  0.47  0.53  0.40  0.25 -0.10 -0.30 -0.53 -0.53 -0.41 -0.15  0.07  0.30  0.46  0.51  0.31  0.16 -0.11 -0.31 -0.44 -0.42
[43] -0.36 -0.13  0.12  0.30  0.37  0.41
 [1]  0.53  0.33 -0.01 -0.31 -0.47 -0.62 -0.50 -0.34  0.03  0.38  0.53  0.54  0.46  0.25 -0.02 -0.28 -0.51 -0.58 -0.50 -0.27 -0.02
[22]  0.28  0.37  0.49  0.41  0.22  0.03 -0.30 -0.48 -0.52 -0.39 -0.23 -0.02  0.24  0.41  0.48  0.42  0.26 -0.01 -0.18 -0.31 -0.43
[43] -0.38 -0.19 -0.01  0.23  0.34  0.40

ccf2(X, Y, 24); ccf2(X, Yw, 24, ylim=c(-.6,.6))

summary(fit <- lm(salmon~time(salmon), na.action=NULL))

Call:
lm(formula = salmon ~ time(salmon), na.action = NULL)

Residuals:
     Min       1Q   Median       3Q      Max 
-1.69187 -0.62453 -0.07024  0.51561  2.34959 

Coefficients:
               Estimate Std. Error t value Pr(>|t|)    
(Intercept)  -503.08947   34.44164  -14.61   <2e-16 ***
time(salmon)    0.25290    0.01713   14.76   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.8814 on 164 degrees of freedom
Multiple R-squared:  0.5706,    Adjusted R-squared:  0.568 
F-statistic: 217.9 on 1 and 164 DF,  p-value: < 2.2e-16

tsplot(salmon, col=4, ylab="USD per KG", main="Salmon Export Price")
abline(fit)

culer = c(rgb(.66,.12,.85), rgb(.12,.66,.85), rgb(.85,.30,.12))
par(mfrow=c(3,1))
tsplot(cmort, main="Cardiovascular Mortality", col=culer[1],

type="o", pch=19, ylab="")

tsplot(tempr, main="Temperature", col=culer[2], type="o", pch=19,

ylab="")

tsplot(part, main="Particulates", col=culer[3], type="o", pch=19,

ylab="")

##-- Figure 3.3 --##
tsplot(cmort, main="", ylab="", ylim=c(20,130), col=culer[1])
lines(tempr, col=culer[2])

lines(part, col=culer[3])
legend("topright", legend=c("Mortality", "Temperature", "Pollution"),

lty=1, lwd=2, col=culer, bg="white")

panel.cor <- function(x, y, ...){
usr <- par("usr"); on.exit(par(usr))
par(usr = c(0, 1, 0, 1))
r <- round(cor(x, y), 2)
text(0.5, 0.5, r, cex = 1.75)
}
pairs(cbind(Mortality=cmort, Temperature=tempr, Particulates=part),

col="dodgerblue3", lower.panel=panel.cor)

par(mfrow = 2:1)
plot(tempr, tempr^2) # collinear
cor(tempr, tempr^2)

temp = tempr - mean(tempr)
plot(temp, temp^2) # not collinear
cor(temp, temp^2)

temp = tempr - mean(tempr) # center temperature
temp2 = temp^2
trend = time(cmort) # time is trend
fit = lm(cmort~ trend + temp + temp2 + part, na.action=NULL)
summary(fit) # regression results

Call:
lm(formula = cmort ~ trend + temp + temp2 + part, na.action = NULL)

Residuals:
     Min       1Q   Median       3Q      Max 
-19.0760  -4.2153  -0.4878   3.7435  29.2448 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)  2.831e+03  1.996e+02   14.19  < 2e-16 ***
trend       -1.396e+00  1.010e-01  -13.82  < 2e-16 ***
temp        -4.725e-01  3.162e-02  -14.94  < 2e-16 ***
temp2        2.259e-02  2.827e-03    7.99 9.26e-15 ***
part         2.554e-01  1.886e-02   13.54  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 6.385 on 503 degrees of freedom
Multiple R-squared:  0.5954,    Adjusted R-squared:  0.5922 
F-statistic:   185 on 4 and 503 DF,  p-value: < 2.2e-16

summary(aov(fit)) # ANOVA table (compare to next line)

             Df Sum Sq Mean Sq F value Pr(>F)    
trend         1  10667   10667  261.62 <2e-16 ***
temp          1   8607    8607  211.09 <2e-16 ***
temp2         1   3429    3429   84.09 <2e-16 ***
part          1   7476    7476  183.36 <2e-16 ***
Residuals   503  20508      41                   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

summary(aov(lm(cmort~cbind(trend, temp, temp2, part)))) # Table 3.1

                                 Df Sum Sq Mean Sq F value Pr(>F)    
cbind(trend, temp, temp2, part)   4  30178    7545     185 <2e-16 ***
Residuals                       503  20508      41                   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

num = length(cmort) # sample size
AIC(fit)/num - log(2*pi) # AIC

[1] 4.721732

BIC(fit)/num - log(2*pi) # BIC

[1] 4.771699

fish = ts.intersect( rec, soiL6=lag(soi,-6) )
summary(fit1 <- lm(rec~ soiL6, data=fish, na.action=NULL))

Call:
lm(formula = rec ~ soiL6, data = fish, na.action = NULL)

Residuals:
    Min      1Q  Median      3Q     Max 
-65.187 -18.234   0.354  16.580  55.790 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   65.790      1.088   60.47   <2e-16 ***
soiL6        -44.283      2.781  -15.92   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 22.5 on 445 degrees of freedom
Multiple R-squared:  0.3629,    Adjusted R-squared:  0.3615 
F-statistic: 253.5 on 1 and 445 DF,  p-value: < 2.2e-16

library(dynlm)
summary(fit2 <- dynlm(rec~ L(soi,6)))

Time series regression with "ts" data:
Start = 1950(7), End = 1987(9)

Call:
dynlm(formula = rec ~ L(soi, 6))

Residuals:
    Min      1Q  Median      3Q     Max 
-65.187 -18.234   0.354  16.580  55.790 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   65.790      1.088   60.47   <2e-16 ***
L(soi, 6)    -44.283      2.781  -15.92   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 22.5 on 445 degrees of freedom
Multiple R-squared:  0.3629,    Adjusted R-squared:  0.3615 
F-statistic: 253.5 on 1 and 445 DF,  p-value: < 2.2e-16

fit = lm(salmon~time(salmon), na.action=NULL) # the regression
par(mfrow=c(2,1)) # plot transformed data
tsplot(resid(fit), main="detrended salmon price")
tsplot(diff(salmon), main="differenced salmon price")

par(mfrow=c(2,1)) # plot their ACFs

acf1(resid(fit), 48, main="detrended salmon price")

 [1]  0.89  0.72  0.56  0.42  0.29  0.20  0.14  0.10  0.08  0.06  0.01 -0.07 -0.20 -0.33 -0.43 -0.50 -0.54 -0.58 -0.56 -0.50 -0.44
[22] -0.35 -0.26 -0.19 -0.17 -0.17 -0.17 -0.17 -0.13 -0.09  0.01  0.14  0.27  0.37  0.45  0.48  0.42  0.34  0.27  0.21  0.17  0.14
[43]  0.14  0.15  0.16  0.14  0.10  0.05

acf1(diff(salmon), 48, main="differenced salmon price")

 [1]  0.26 -0.02 -0.08 -0.10 -0.19 -0.11 -0.10 -0.06  0.02  0.07  0.15  0.24  0.00 -0.15 -0.14 -0.07 -0.04 -0.23 -0.15 -0.05 -0.07
[22] -0.06  0.13  0.21  0.11 -0.01  0.00 -0.16 -0.09 -0.23 -0.10  0.04  0.12  0.09  0.23  0.34  0.08 -0.06 -0.06 -0.08 -0.07 -0.09
[43] -0.06  0.06  0.10  0.06  0.04  0.12

par(mfrow=c(2,1))
tsplot(diff(gtemp_land), col=4, main="differenced global temperature")
mean(diff(gtemp_land)) # drift since 1880

[1] 0.01595376

acf1(diff(gtemp_land))

 [1] -0.51  0.04  0.00 -0.05  0.00  0.09 -0.06 -0.06  0.13 -0.05 -0.13  0.24 -0.15  0.01  0.07 -0.05 -0.02  0.03  0.05 -0.08  0.03
[22]  0.00 -0.08  0.14

mean(window(diff(gtemp_land), start=1980)) # drift since 1980

[1] 0.04909091

layout(matrix(1:4,2), widths=c(2.5,1))
par(mgp=c(1.6,.6,0), mar=c(2,2,.5,0)+.5)
tsplot(varve, main="", ylab="", col=4, margin=0)
mtext("varve", side=3, line=.5, cex=1.2, font=2, adj=0)

tsplot(log(varve), main="", ylab="", col=4, margin=0)
mtext("log(varve)", side=3, line=.5, cex=1.2, font=2, adj=0)

qqnorm(varve, main="", col=4); qqline(varve, col=2, lwd=2)
qqnorm(log(varve), main="", col=4); qqline(log(varve), col=2, lwd=2)

lag1.plot(soi, 12, col="dodgerblue3") # Figure 3.10

lag2.plot(soi, rec, 8, col="dodgerblue3") # Figure 3.11

library(zoo) # zoo allows easy use of the variable names
dummy = ifelse(soi<0, 0, 1)
fish = as.zoo(ts.intersect(rec, soiL6=lag(soi,-6), dL6=lag(dummy,-6)))
summary(fit <- lm(rec~ soiL6*dL6, data=fish, na.action=NULL))

Call:
lm(formula = rec ~ soiL6 * dL6, data = fish, na.action = NULL)

Residuals:
    Min      1Q  Median      3Q     Max 
-63.291 -15.821   2.224  15.791  61.788 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   74.479      2.865  25.998  < 2e-16 ***
soiL6        -15.358      7.401  -2.075   0.0386 *  
dL6           -1.139      3.711  -0.307   0.7590    
soiL6:dL6    -51.244      9.523  -5.381  1.2e-07 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 21.84 on 443 degrees of freedom
Multiple R-squared:  0.4024,    Adjusted R-squared:  0.3984 
F-statistic: 99.43 on 3 and 443 DF,  p-value: < 2.2e-16

plot(fish$soiL6, fish$rec, panel.first=Grid(), col="dodgerblue3")
points(fish$soiL6, fitted(fit), pch=3, col=6)

lines(lowess(fish$soiL6, fish$rec), col=4, lwd=2)

tsplot(resid(fit)) # not shown, but looks like Figure 3.5

acf1(resid(fit)) # and obviously not noise

 [1]  0.69  0.62  0.49  0.37  0.24  0.15  0.08  0.00 -0.03 -0.10 -0.13 -0.16 -0.17 -0.23 -0.24 -0.23 -0.23 -0.22 -0.17 -0.09 -0.05
[22]  0.01  0.05  0.06  0.09  0.07  0.10  0.06  0.02 -0.02 -0.02 -0.02

set.seed(90210) # so you can reproduce these results
x = 2*cos(2*pi*1:500/50 + .6*pi) + rnorm(500,0,5)
z1 = cos(2*pi*1:500/50)
z2 = sin(2*pi*1:500/50)
summary(fit <- lm(x~ 0 + z1 + z2)) # zero to exclude the intercept

Call:
lm(formula = x ~ 0 + z1 + z2)

Residuals:
     Min       1Q   Median       3Q      Max 
-14.8584  -3.8525  -0.3186   3.3487  15.5440 

Coefficients:
   Estimate Std. Error t value Pr(>|t|)    
z1  -0.7442     0.3274  -2.273   0.0235 *  
z2  -1.9949     0.3274  -6.093 2.23e-09 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 5.177 on 498 degrees of freedom
Multiple R-squared:  0.07827,   Adjusted R-squared:  0.07456 
F-statistic: 21.14 on 2 and 498 DF,  p-value: 1.538e-09

par(mfrow=c(2,1))
tsplot(x, col=4)
tsplot(x, ylab=expression(hat(x)), col=rgb(0,0,1,.5))

lines(fitted(fit), col=2, lwd=2)

w = c(.5, rep(1,11), .5)/12
soif = filter(soi, sides=2, filter=w)
tsplot(soi, col=rgb(.5, .6, .85, .9), ylim=c(-1, 1.15))
lines(soif, lwd=2, col=4)

# insert
par(fig = c(.65, 1, .75, 1), new = TRUE)
w1 = c(rep(0,20), w, rep(0,20))
plot(w1, type="l", ylim = c(-.02,.1), xaxt="n", yaxt="n", ann=FALSE)

tsplot(soi, col=rgb(0.5, 0.6, 0.85, .9), ylim=c(-1, 1.15))
lines(ksmooth(time(soi), soi, "normal", bandwidth=1), lwd=2, col=4)

# insert
par(fig = c(.65, 1, .75, 1), new = TRUE)
curve(dnorm(x), -3, 3, xaxt="n", yaxt="n", ann=FALSE, col=4)

SOI = ts(soi, freq=1)
tsplot(SOI) # the time scale matters (not shown)
lines(ksmooth(time(SOI), SOI, "normal", bandwidth=12), lwd=2, col=4)

tsplot(soi, col=rgb(0.5, 0.6, 0.85, .9))
lines(lowess(soi, f=.05), lwd=2, col=4) # El Niño cycle

# lines(lowess(soi), lty=2, lwd=2, col=2) # trend (with default span)
##-- trend with CIs using loess --##
lo = predict(loess(soi~ time(soi)), se=TRUE)
trnd = ts(lo$fit, start=1950, freq=12) # put back ts attributes
lines(trnd, col=6, lwd=2)
L = trnd - qt(0.975, lo$df)*lo$se
U = trnd + qt(0.975, lo$df)*lo$se
xx = c(time(soi), rev(time(soi)))
yy = c(L, rev(U))
polygon(xx, yy, border=8, col=gray(.6, alpha=.4) )

plot(tempr, cmort, xlab="Temperature", ylab="Mortality",
col="dodgerblue3", panel.first=Grid())
lines(lowess(tempr,cmort), col=4, lwd=2)

x = window(hor, start=2002)
plot(decompose(x)) # not shown

plot(stl(x, s.window="per")) # seasons are periodic - not shown

plot(stl(x, s.window=15))

culer = c("cyan4", 4, 2, 6)
par(mfrow = c(4,1), cex.main=1)
x = window(hor, start=2002)
out = stl(x, s.window=15)$time.series
tsplot(x, main="Hawaiian Occupancy Rate", ylab="% rooms", col=gray(.7))
text(x, labels=1:4, col=culer, cex=1.25)

tsplot(out[,1], main="Seasonal", ylab="% rooms",col=gray(.7))
text(out[,1], labels=1:4, col=culer, cex=1.25)

tsplot(out[,2], main="Trend", ylab="% rooms", col=gray(.7))
text(out[,2], labels=1:4, col=culer, cex=1.25)

tsplot(out[,3], main="Noise", ylab="% rooms", col=gray(.7))
text(out[,3], labels=1:4, col=culer, cex=1.25)

par(mfrow=c(2,1))
tsplot(arima.sim(list(order=c(1,0,0), ar=.9), n=100), ylab="x", col=4,

main=expression(AR(1)~~~phi==+.9))

tsplot(arima.sim(list(order=c(1,0,0),ar=-.9), n=100), ylab="x", col=4,

main=expression(AR(1)~~~phi==-.9))

psi = ARMAtoMA(ar = c(1.5, -.75), ma = 0, 50)
par(mfrow=c(2,1), mar=c(2,2.5,1,0)+.5, mgp=c(1.5,.6,0), cex.main=1.1)
plot(psi, xaxp=c(0,144,12), type="n", col=4,
ylab=expression(psi-weights),
main=expression(AR(2)~~~phi[1]==1.5~~~phi[2]==-.75))
abline(v=seq(0,48,by=12), h=seq(-.5,1.5,.5), col=gray(.9))

lines(psi, type="o", col=4)

set.seed(8675309)
simulation = arima.sim(list(order=c(2,0,0), ar=c(1.5,-.75)), n=144)

plot(simulation, xaxp=c(0,144,12), type="n", ylab=expression(X[~t]))
abline(v=seq(0,144,by=12), h=c(-5,0,5), col=gray(.9))

lines(simulation, col=4)

par(mfrow = c(2,1))
tsplot(arima.sim(list(order=c(0,0,1), ma=.9), n=100), col=4,
ylab="x", main=expression(MA(1)~~~theta==+.5))
tsplot(arima.sim(list(order=c(0,0,1), ma=-.9), n=100), col=4,
ylab="x", main=expression(MA(1)~~~theta==-.5))

set.seed(8675309) # Jenny, I got your number
x = rnorm(150, mean=5) # generate iid N(5,1)s
arima(x, order=c(1,0,1)) # estimation

Call:
arima(x = x, order = c(1, 0, 1))

Coefficients:
          ar1     ma1  intercept
      -0.9595  0.9527     5.0462
s.e.   0.1688  0.1750     0.0727

sigma^2 estimated as 0.7986:  log likelihood = -195.98,  aic = 399.96

AR = c(1, -.3, -.4) # original AR coefs on the left
polyroot(AR)

[1]  1.25-0i -2.00+0i

MA = c(1, .5) # original MA coefs on the right
polyroot(MA)

[1] -2+0i

ACF = ARMAacf(ar=c(1.5,-.75), ma=0, 50)
plot(ACF, type="h", xlab="lag", panel.first=Grid())
abline(h=0)

ACF = ARMAacf(ar=c(1.5,-.75), ma=0, 24)[-1]
PACF = ARMAacf(ar=c(1.5,-.75), ma=0, 24, pacf=TRUE)
par(mfrow=1:2)
tsplot(ACF, type="h", xlab="lag", ylim=c(-.8,1))
abline(h=0)

tsplot(PACF, type="h", xlab="lag", ylim=c(-.8,1))
abline(h=0)

acf2(rec, 48) # will produce values and a graphic

     [,1]  [,2]  [,3]  [,4] [,5]  [,6]  [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21]
ACF  0.92  0.78  0.63  0.48 0.36  0.26  0.18 0.13 0.09  0.07  0.06  0.02 -0.04 -0.12 -0.19 -0.24 -0.27 -0.27 -0.24 -0.19 -0.11
PACF 0.92 -0.44 -0.05 -0.02 0.07 -0.03 -0.03 0.04 0.05 -0.02 -0.05 -0.14 -0.15 -0.05  0.05  0.01  0.01  0.02  0.09  0.11  0.03
     [,22] [,23] [,24] [,25] [,26] [,27] [,28] [,29] [,30] [,31] [,32] [,33] [,34] [,35] [,36] [,37] [,38] [,39] [,40] [,41]
ACF  -0.03  0.03  0.06  0.06  0.02 -0.02 -0.06 -0.09 -0.12 -0.13 -0.11 -0.05  0.02  0.08  0.12  0.10  0.06  0.01 -0.02 -0.03
PACF -0.03 -0.01 -0.07 -0.12 -0.03  0.05 -0.08 -0.04 -0.03  0.06  0.05  0.15  0.09 -0.04 -0.10 -0.09 -0.02  0.05  0.08 -0.02
     [,42] [,43] [,44] [,45] [,46] [,47] [,48]
ACF  -0.03 -0.02  0.01  0.06  0.12  0.17  0.20
PACF -0.01 -0.02  0.05  0.01  0.05  0.08 -0.04

(regr = ar.ols(rec, order=2, demean=FALSE, intercept=TRUE))

Call:
ar.ols(x = rec, order.max = 2, demean = FALSE, intercept = TRUE)

Coefficients:
      1        2  
 1.3541  -0.4632  

Intercept: 6.737 (1.111) 

Order selected 2  sigma^2 estimated as  89.72

regr$asy.se.coef # standard errors of the estimates

$x.mean
[1] 1.110599

$ar
[1] 0.04178901 0.04187942

set.seed(2)
ma1 = arima.sim(list(order = c(0,0,1), ma = 0.9), n = 50)
acf1(ma1, plot=FALSE)[1]

[1] 0.5066599

tsplot(diff(log(varve)), col=4, ylab=expression(nabla~log~X[~t]),

main="Transformed Glacial Varves")

acf2(diff(log(varve)))

     [,1]  [,2]  [,3]  [,4]  [,5]  [,6]  [,7]  [,8]  [,9] [,10] [,11] [,12] [,13] [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21]
ACF  -0.4 -0.04 -0.06  0.01  0.00  0.04 -0.04  0.04  0.01 -0.05  0.06 -0.06 -0.04  0.08 -0.02  0.01  0.00  0.03 -0.05 -0.06  0.07
PACF -0.4 -0.24 -0.23 -0.18 -0.15 -0.08 -0.11 -0.05 -0.01 -0.07  0.02 -0.05 -0.12 -0.03 -0.05 -0.04 -0.03  0.03 -0.03 -0.13 -0.04
     [,22] [,23] [,24] [,25] [,26] [,27] [,28] [,29] [,30] [,31] [,32] [,33] [,34] [,35] [,36]
ACF   0.04 -0.06  0.05 -0.01 -0.04  0.05 -0.05  0.03 -0.02  0.00  0.06 -0.05 -0.03  0.04 -0.05
PACF  0.01 -0.06  0.01  0.02 -0.04  0.03 -0.02  0.00 -0.03 -0.02  0.04 -0.02 -0.02  0.01 -0.07

x = diff(log(varve)) # data
r = acf1(x, 1, plot=FALSE) # acf(1)
c(0) -> w -> z -> Sc -> Sz -> Szw -> para # initialize
num = length(x) # = 633
## Estimation
para[1] = (1-sqrt(1-4*(r^2)))/(2*r) # MME
niter = 12
for (j in 1:niter){
for (i in 2:num){ w[i] = x[i] - para[j]*w[i-1]
z[i] = w[i-1] - para[j]*z[i-1]

}
Sc[j] = sum(w^2)
Sz[j] = sum(z^2)
Szw[j] = sum(z*w)
para[j+1] = para[j] + Szw[j]/Sz[j]
}
## Results
cbind(iteration=1:niter-1, thetahat=para[1:niter], Sc, Sz)

      iteration   thetahat       Sc       Sz
 [1,]         0 -0.4946886 158.7393 171.2397
 [2,]         1 -0.6681759 150.7468 235.2660
 [3,]         2 -0.7333163 149.2644 300.5618
 [4,]         3 -0.7563893 149.0309 336.8234
 [5,]         4 -0.7655639 148.9897 354.1729
 [6,]         5 -0.7694700 148.9818 362.1665
 [7,]         6 -0.7711856 148.9803 365.8014
 [8,]         7 -0.7719497 148.9800 367.4456
 [9,]         8 -0.7722921 148.9799 368.1875
[10,]         9 -0.7724460 148.9799 368.5220
[11,]        10 -0.7725153 148.9799 368.6727
[12,]        11 -0.7725465 148.9799 368.7406

## Plot conditional SS
c(0) -> w -> cSS
th = -seq(.3, .94, .01)
for (p in 1:length(th)){
for (i in 2:num){ w[i] = x[i] - th[p]*w[i-1]
}
cSS[p] = sum(w^2)
}
tsplot(th, cSS, ylab=expression(S[c](theta)), xlab=expression(theta))
abline(v=para[1:12], lty=2, col=4) # add previous results to plot

points(para[1:12], Sc[1:12], pch=16, col=4)

sarima(diff(log(varve)), p=0, d=0, q=1, no.constant=TRUE)

initial  value -0.551778 
iter   2 value -0.671626
iter   3 value -0.705973
iter   4 value -0.707314
iter   5 value -0.722372
iter   6 value -0.722738
iter   7 value -0.723187
iter   8 value -0.723194
iter   9 value -0.723195
iter   9 value -0.723195
iter   9 value -0.723195
final  value -0.723195 
converged
initial  value -0.722700 
iter   2 value -0.722702
iter   3 value -0.722702
iter   3 value -0.722702
iter   3 value -0.722702
final  value -0.722702 
converged
<><><><><><><><><><><><><><>
 
Coefficients: 
    Estimate     SE  t.value p.value
ma1  -0.7705 0.0341 -22.6161       0

sigma^2 estimated as 0.2353156 on 632 degrees of freedom 
 
AIC = 1.398791  AICc = 1.398802  BIC = 1.412853 
 

sarima(rec, p=2, d=0, q=0)

initial  value 3.332380 
iter   2 value 3.251366
iter   3 value 2.564654
iter   4 value 2.430141
iter   5 value 2.258212
iter   6 value 2.253343
iter   7 value 2.248346
iter   8 value 2.248345
iter   9 value 2.248345
iter  10 value 2.248341
iter  11 value 2.248332
iter  12 value 2.248331
iter  13 value 2.248330
iter  13 value 2.248330
iter  13 value 2.248330
final  value 2.248330 
converged
initial  value 2.248862 
iter   2 value 2.248857
iter   3 value 2.248855
iter   4 value 2.248855
iter   5 value 2.248854
iter   6 value 2.248854
iter   7 value 2.248854
iter   8 value 2.248854
iter   9 value 2.248853
iter  10 value 2.248853
iter  10 value 2.248853
iter  10 value 2.248853
final  value 2.248853 
converged
<><><><><><><><><><><><><><>
 
Coefficients: 
      Estimate     SE  t.value p.value
ar1     1.3512 0.0416  32.4933       0
ar2    -0.4612 0.0417 -11.0687       0
xmean  61.8585 4.0039  15.4494       0

sigma^2 estimated as 89.33436 on 450 degrees of freedom 
 
AIC = 7.353244  AICc = 7.353362  BIC = 7.389587 
 

sarima.for(rec, n.ahead=24, p=2, d=0, q=0)

$pred
          Jan      Feb      Mar      Apr      May      Jun      Jul      Aug      Sep      Oct      Nov      Dec
1987                                                                                  20.36547 26.08036 32.65148
1988 38.89474 44.30006 48.72437 52.20958 54.87831 56.87693 58.34666 59.41079 60.17081 60.70697 61.08091 61.33890
1989 61.51504 61.63405 61.71363 61.76626 61.80068 61.82291 61.83707 61.84596 61.85143                           

$se
           Jan       Feb       Mar       Apr       May       Jun       Jul       Aug       Sep       Oct       Nov       Dec
1987                                                                                            9.451686 15.888378 20.464325
1988 23.492457 25.393693 26.537088 27.199368 27.570234 27.771616 27.877923 27.932597 27.960042 27.973508 27.979974 27.983014
1989 27.984413 27.985045 27.985323 27.985444 27.985495 27.985515 27.985524 27.985527 27.985528                              

abline(h=61.8585, col=4) # display estimated mean

sarima(diff(log(varve)), p=0, d=0, q=1, no.constant=TRUE)

initial  value -0.551778 
iter   2 value -0.671626
iter   3 value -0.705973
iter   4 value -0.707314
iter   5 value -0.722372
iter   6 value -0.722738
iter   7 value -0.723187
iter   8 value -0.723194
iter   9 value -0.723195
iter   9 value -0.723195
iter   9 value -0.723195
final  value -0.723195 
converged
initial  value -0.722700 
iter   2 value -0.722702
iter   3 value -0.722702
iter   3 value -0.722702
iter   3 value -0.722702
final  value -0.722702 
converged
<><><><><><><><><><><><><><>
 
Coefficients: 
    Estimate     SE  t.value p.value
ma1  -0.7705 0.0341 -22.6161       0

sigma^2 estimated as 0.2353156 on 632 degrees of freedom 
 
AIC = 1.398791  AICc = 1.398802  BIC = 1.412853 
 

sarima(log(varve), p=0, d=1, q=1, no.constant=TRUE)

initial  value -0.551778 
iter   2 value -0.671626
iter   3 value -0.705973
iter   4 value -0.707314
iter   5 value -0.722372
iter   6 value -0.722738
iter   7 value -0.723187
iter   8 value -0.723194
iter   9 value -0.723195
iter   9 value -0.723195
iter   9 value -0.723195
final  value -0.723195 
converged
initial  value -0.722700 
iter   2 value -0.722702
iter   3 value -0.722702
iter   3 value -0.722702
iter   3 value -0.722702
final  value -0.722702 
converged
<><><><><><><><><><><><><><>
 
Coefficients: 
    Estimate     SE  t.value p.value
ma1  -0.7705 0.0341 -22.6161       0

sigma^2 estimated as 0.2353156 on 632 degrees of freedom 
 
AIC = 1.398792  AICc = 1.398802  BIC = 1.412853 
 

ARMAtoMA(ar=1, ma=0, 20) # ψ-weights

 [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1

round( ARMAtoMA(ar=c(1.9,-.9), ma=0, 60), 1 )

 [1]  1.9  2.7  3.4  4.1  4.7  5.2  5.7  6.1  6.5  6.9  7.2  7.5  7.7  7.9  8.1  8.3  8.5  8.6  8.8  8.9  9.0  9.1  9.2  9.3  9.4
[26]  9.4  9.5  9.5  9.6  9.6  9.7  9.7  9.7  9.7  9.8  9.8  9.8  9.8  9.9  9.9  9.9  9.9  9.9  9.9  9.9  9.9  9.9  9.9  9.9 10.0
[51] 10.0 10.0 10.0 10.0 10.0 10.0 10.0 10.0 10.0 10.0

set.seed(1998)
x <- ts(arima.sim(list(order = c(1,1,0), ar=.9), n=150)[-1])
y <- window(x, start=1, end=100)
sarima.for(y, n.ahead = 50, p = 1, d = 1, q = 0, plot.all=TRUE)

$pred
Time Series:
Start = 101 
End = 150 
Frequency = 1 
 [1] 108.8050 107.6806 106.8693 106.3241 106.0053 105.8791 105.9165 106.0932 106.3883 106.7841 107.2656 107.8199 108.4361 109.1050
[15] 109.8187 110.5705 111.3547 112.1665 113.0017 113.8568 114.7289 115.6154 116.5141 117.4233 118.3413 119.2669 120.1989 121.1363
[29] 122.0784 123.0244 123.9738 124.9260 125.8807 126.8374 127.7959 128.7558 129.7171 130.6794 131.6426 132.6066 133.5713 134.5365
[43] 135.5022 136.4684 137.4348 138.4016 139.3686 140.3358 141.3033 142.2708

$se
Time Series:
Start = 101 
End = 150 
Frequency = 1 
 [1]  0.9570915  2.0131097  3.1812332  4.4084707  5.6617584  6.9197429  8.1683744  9.3983950 10.6037992 11.7808348 12.9273310
[12] 14.0422326 15.1252705 16.1767257 17.1972588 18.1877853 19.1493851 20.0832363 20.9905665 21.8726189 22.7306265 23.5657957
[23] 24.3792934 25.1722400 25.9457044 26.7007014 27.4381914 28.1590801 28.8642201 29.5544126 30.2304096 30.8929162 31.5425932
[34] 32.1800597 32.8058954 33.4206433 34.0248122 34.6188788 35.2032899 35.7784645 36.3447959 36.9026532 37.4523833 37.9943122
[45] 38.5287468 39.0559758 39.5762716 40.0898907 40.5970755 41.0980549

text(85, 205, "PAST"); text(115, 205, "FUTURE")

abline(v=100, lty=2, col=4)
lines(x)

set.seed(666)
x = arima.sim(list(order = c(0,1,1), ma = -0.8), n = 100)
(x.ima = HoltWinters(x, beta=FALSE, gamma=FALSE)) # α below is 1 − λ

Holt-Winters exponential smoothing without trend and without seasonal component.

Call:
HoltWinters(x = x, beta = FALSE, gamma = FALSE)

Smoothing parameters:
 alpha: 0.1663072
 beta : FALSE
 gamma: FALSE

Coefficients:
       [,1]
a -2.241533

plot(x.ima, main="EWMA")

##-- Figure 5.3 --##
layout(1:2, heights=2:1)
tsplot(gnp, col=4)
acf1(gnp, main="")

 [1] 0.99 0.97 0.96 0.94 0.93 0.91 0.90 0.88 0.87 0.85 0.83 0.82 0.80 0.79 0.77 0.76 0.74 0.73 0.72 0.70 0.69 0.68 0.66 0.65 0.64

##-- Figure 5.4 --##
tsplot(diff(log(gnp)), ylab="GNP Growth Rate", col=4)

sarima(diff(log(gnp)), 0,0,2) # MA(2) on growth rate

initial  value -4.591629 
iter   2 value -4.661095
iter   3 value -4.662220
iter   4 value -4.662243
iter   5 value -4.662243
iter   6 value -4.662243
iter   6 value -4.662243
iter   6 value -4.662243
final  value -4.662243 
converged
initial  value -4.662022 
iter   2 value -4.662023
iter   2 value -4.662023
iter   2 value -4.662023
final  value -4.662023 
converged
<><><><><><><><><><><><><><>
 
Coefficients: 
      Estimate     SE t.value p.value
ma1     0.3028 0.0654  4.6272  0.0000
ma2     0.2035 0.0644  3.1594  0.0018
xmean   0.0083 0.0010  8.7178  0.0000

sigma^2 estimated as 8.919178e-05 on 219 degrees of freedom 
 
AIC = -6.450133  AICc = -6.449637  BIC = -6.388823 
 

sarima(diff(log(gnp)), 1,0,0)

initial  value -4.589567 
iter   2 value -4.654150
iter   3 value -4.654150
iter   4 value -4.654151
iter   4 value -4.654151
iter   4 value -4.654151
final  value -4.654151 
converged
initial  value -4.655919 
iter   2 value -4.655921
iter   3 value -4.655922
iter   4 value -4.655922
iter   5 value -4.655922
iter   5 value -4.655922
iter   5 value -4.655922
final  value -4.655922 
converged
<><><><><><><><><><><><><><>
 
Coefficients: 
      Estimate     SE t.value p.value
ar1     0.3467 0.0627  5.5255       0
xmean   0.0083 0.0010  8.5398       0

sigma^2 estimated as 9.029569e-05 on 220 degrees of freedom 
 
AIC = -6.44694  AICc = -6.446693  BIC = -6.400958 
 

round( ARMAtoMA(ar=.35, ma=0, 10), 3) # print psi-weights

 [1] 0.350 0.122 0.043 0.015 0.005 0.002 0.001 0.000 0.000 0.000

sarima(diff(log(gnp)), 0, 0, 2) # MA(2) fit with diagnostics

initial  value -4.591629 
iter   2 value -4.661095
iter   3 value -4.662220
iter   4 value -4.662243
iter   5 value -4.662243
iter   6 value -4.662243
iter   6 value -4.662243
iter   6 value -4.662243
final  value -4.662243 
converged
initial  value -4.662022 
iter   2 value -4.662023
iter   2 value -4.662023
iter   2 value -4.662023
final  value -4.662023 
converged
<><><><><><><><><><><><><><>
 
Coefficients: 
      Estimate     SE t.value p.value
ma1     0.3028 0.0654  4.6272  0.0000
ma2     0.2035 0.0644  3.1594  0.0018
xmean   0.0083 0.0010  8.7178  0.0000

sigma^2 estimated as 8.919178e-05 on 219 degrees of freedom 
 
AIC = -6.450133  AICc = -6.449637  BIC = -6.388823 
 

sarima(diff(log(gnp)), 0, 0, 3) # try an MA(2+1) fit (not shown)

initial  value -4.591629 
iter   2 value -4.660639
iter   3 value -4.665404
iter   4 value -4.665858
iter   5 value -4.665953
iter   6 value -4.665955
iter   7 value -4.665956
iter   8 value -4.665956
iter   8 value -4.665956
iter   8 value -4.665956
final  value -4.665956 
converged
initial  value -4.665676 
iter   2 value -4.665678
iter   3 value -4.665678
iter   3 value -4.665678
iter   3 value -4.665678
final  value -4.665678 
converged
<><><><><><><><><><><><><><>
 
Coefficients: 
      Estimate     SE t.value p.value
ma1     0.3208 0.0662  4.8430  0.0000
ma2     0.2478 0.0718  3.4512  0.0007
ma3     0.0909 0.0701  1.2962  0.1963
xmean   0.0083 0.0010  7.9444  0.0000

sigma^2 estimated as 8.852732e-05 on 218 degrees of freedom 
 
AIC = -6.448434  AICc = -6.447604  BIC = -6.371797 
 

sarima(diff(log(gnp)), 2, 0, 0) # try an AR(1+1) fit (not shown)

initial  value -4.588072 
iter   2 value -4.647096
iter   3 value -4.655676
iter   4 value -4.655896
iter   5 value -4.655903
iter   6 value -4.655903
iter   7 value -4.655904
iter   7 value -4.655904
iter   7 value -4.655904
final  value -4.655904 
converged
initial  value -4.659296 
iter   2 value -4.659300
iter   3 value -4.659306
iter   4 value -4.659307
iter   5 value -4.659307
iter   5 value -4.659307
iter   5 value -4.659307
final  value -4.659307 
converged
<><><><><><><><><><><><><><>
 
Coefficients: 
      Estimate     SE t.value p.value
ar1     0.3180 0.0667  4.7660  0.0000
ar2     0.0820 0.0668  1.2282  0.2207
xmean   0.0083 0.0011  7.8742  0.0000

sigma^2 estimated as 8.968109e-05 on 219 degrees of freedom 
 
AIC = -6.444701  AICc = -6.444205  BIC = -6.383392 
 

sarima(log(varve), 0, 1, 1, no.constant=TRUE) # ARIMA(0,1,1)

initial  value -0.551778 
iter   2 value -0.671626
iter   3 value -0.705973
iter   4 value -0.707314
iter   5 value -0.722372
iter   6 value -0.722738
iter   7 value -0.723187
iter   8 value -0.723194
iter   9 value -0.723195
iter   9 value -0.723195
iter   9 value -0.723195
final  value -0.723195 
converged
initial  value -0.722700 
iter   2 value -0.722702
iter   3 value -0.722702
iter   3 value -0.722702
iter   3 value -0.722702
final  value -0.722702 
converged
<><><><><><><><><><><><><><>
 
Coefficients: 
    Estimate     SE  t.value p.value
ma1  -0.7705 0.0341 -22.6161       0

sigma^2 estimated as 0.2353156 on 632 degrees of freedom 
 
AIC = 1.398792  AICc = 1.398802  BIC = 1.412853 
 

sarima(log(varve), 1, 1, 1, no.constant=TRUE) # ARIMA(1,1,1)

initial  value -0.550992 
iter   2 value -0.648952
iter   3 value -0.676952
iter   4 value -0.699136
iter   5 value -0.724481
iter   6 value -0.726964
iter   7 value -0.734257
iter   8 value -0.735999
iter   9 value -0.737045
iter  10 value -0.737381
iter  11 value -0.737469
iter  12 value -0.737473
iter  13 value -0.737473
iter  14 value -0.737473
iter  14 value -0.737473
iter  14 value -0.737473
final  value -0.737473 
converged
initial  value -0.737355 
iter   2 value -0.737361
iter   3 value -0.737362
iter   4 value -0.737363
iter   5 value -0.737363
iter   5 value -0.737363
iter   5 value -0.737363
final  value -0.737363 
converged
<><><><><><><><><><><><><><>
 
Coefficients: 
    Estimate     SE  t.value p.value
ar1   0.2330 0.0518   4.4994       0
ma1  -0.8858 0.0292 -30.3861       0

sigma^2 estimated as 0.2284339 on 631 degrees of freedom 
 
AIC = 1.37263  AICc = 1.372661  BIC = 1.393723 
 

uspop = c(75.995, 91.972, 105.711, 123.203, 131.669,150.697,
179.323, 203.212, 226.505, 249.633, 281.422, 308.745)
uspop = ts(uspop, start=1900, freq=.1)
t = time(uspop) - 1955
reg = lm( uspop~ t+I(t^2)+I(t^3)+I(t^4)+I(t^5)+I(t^6)+I(t^7)+I(t^8) )

b = as.vector(reg$coef)
g = function(t){ b[1] + b[2]*(t-1955) + b[3]*(t-1955)^2 +
b[4]*(t-1955)^3 + b[5]*(t-1955)^4 + b[6]*(t-1955)^5 +
b[7]*(t-1955)^6 + b[8]*(t-1955)^7 + b[9]*(t-1955)^8

}
par(mar=c(2,2.5,.5,0)+.5, mgp=c(1.6,.6,0))
curve(g, 1900, 2024, ylab="Population", xlab="Year", main="U.S.
Population by Official Census", panel.first=Grid(),
cex.main=1, font.main=1, col=4)
abline(v=seq(1910,2020,by=20), lty=1, col=gray(.9))

points(time(uspop), uspop, pch=21, bg=rainbow(12), cex=1.25)
mtext(expression(""%*% 10^6), side=2, line=1.5, adj=.95)

axis(1, seq(1910,2020,by=20), labels=TRUE)

sarima(diff(log(gnp)), 1, 0, 0) # AR(1)

initial  value -4.589567 
iter   2 value -4.654150
iter   3 value -4.654150
iter   4 value -4.654151
iter   4 value -4.654151
iter   4 value -4.654151
final  value -4.654151 
converged
initial  value -4.655919 
iter   2 value -4.655921
iter   3 value -4.655922
iter   4 value -4.655922
iter   5 value -4.655922
iter   5 value -4.655922
iter   5 value -4.655922
final  value -4.655922 
converged
<><><><><><><><><><><><><><>
 
Coefficients: 
      Estimate     SE t.value p.value
ar1     0.3467 0.0627  5.5255       0
xmean   0.0083 0.0010  8.5398       0

sigma^2 estimated as 9.029569e-05 on 220 degrees of freedom 
 
AIC = -6.44694  AICc = -6.446693  BIC = -6.400958 
 

sarima(diff(log(gnp)), 0, 0, 2) # MA(2)

initial  value -4.591629 
iter   2 value -4.661095
iter   3 value -4.662220
iter   4 value -4.662243
iter   5 value -4.662243
iter   6 value -4.662243
iter   6 value -4.662243
iter   6 value -4.662243
final  value -4.662243 
converged
initial  value -4.662022 
iter   2 value -4.662023
iter   2 value -4.662023
iter   2 value -4.662023
final  value -4.662023 
converged
<><><><><><><><><><><><><><>
 
Coefficients: 
      Estimate     SE t.value p.value
ma1     0.3028 0.0654  4.6272  0.0000
ma2     0.2035 0.0644  3.1594  0.0018
xmean   0.0083 0.0010  8.7178  0.0000

sigma^2 estimated as 8.919178e-05 on 219 degrees of freedom 
 
AIC = -6.450133  AICc = -6.449637  BIC = -6.388823 
 

set.seed(666)
phi = c(rep(0,11),.9)
sAR = ts(arima.sim(list(order=c(12,0,0), ar=phi), n=37), freq=12) + 50
layout(matrix(c(1,2, 1,3), nc=2), heights=c(1.5,1))
par(mar=c(2.5,2.5,2,1), mgp=c(1.6,.6,0))
plot(sAR, xaxt="n", col=gray(.6), main="seasonal AR(1)", xlab="YEAR",

type="c", ylim=c(45,54))
abline(v=1:4, lty=2, col=gray(.6))

axis(1,1:4); box()
abline(h=seq(46,54,by=2), col=gray(.9))

Months = c("J","F","M","A","M","J","J","A","S","O","N","D")
points(sAR, pch=Months, cex=1.35, font=4, col=1:4)
ACF = ARMAacf(ar=phi, ma=0, 100)[-1]
PACF = ARMAacf(ar=phi, ma=0, 100, pacf=TRUE)
LAG = 1:100/12
plot(LAG, ACF, type="h", xlab="LAG", ylim=c(-.1,1), axes=FALSE)

segments(0,0,0,1)

axis(1, seq(0,8,by=1)); axis(2); box(); abline(h=0)

plot(LAG, PACF, type="h", xlab="LAG", ylim=c(-.1,1), axes=FALSE)
axis(1, seq(0,8,by=1)); axis(2); box(); abline(h=0)

##-- Figure 5.10 --##
phi = c(rep(0,11),.8)
ACF = ARMAacf(ar=phi, ma=-.5, 50)[-1]
PACF = ARMAacf(ar=phi, ma=-.5, 50, pacf=TRUE)
LAG = 1:50/12
par(mfrow=c(1,2))
plot(LAG, ACF, type="h", ylim=c(-.4,.8), panel.first=Grid())
abline(h=0)

plot(LAG, PACF, type="h", ylim=c(-.4,.8), panel.first=Grid())
abline(h=0)

##-- birth series --##
tsplot(birth) # monthly number of births in US

acf2( diff(birth) ) # P/ACF of the differenced birth rate

      [,1] [,2]  [,3]  [,4]  [,5]  [,6]  [,7]  [,8]  [,9] [,10] [,11] [,12] [,13] [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21]
ACF  -0.32 0.16 -0.08 -0.19  0.09 -0.28  0.06 -0.19 -0.05  0.17 -0.26  0.82 -0.28  0.17 -0.07 -0.18  0.08 -0.28  0.07 -0.18 -0.05
PACF -0.32 0.06 -0.01 -0.25 -0.03 -0.26 -0.17 -0.29 -0.35 -0.16 -0.59  0.57  0.13  0.11  0.13  0.09  0.00  0.00  0.05  0.04 -0.07
     [,22] [,23] [,24] [,25] [,26] [,27] [,28] [,29] [,30] [,31] [,32] [,33] [,34] [,35] [,36] [,37] [,38] [,39] [,40] [,41]
ACF   0.16 -0.24  0.78 -0.27  0.19 -0.08 -0.17  0.07 -0.29  0.07 -0.15 -0.04  0.14 -0.24  0.75 -0.23  0.16 -0.08 -0.15  0.05
PACF -0.10 -0.20  0.19  0.01  0.05  0.07  0.07 -0.02 -0.06 -0.02  0.09  0.03 -0.06 -0.16  0.03  0.08 -0.10 -0.03  0.07 -0.04
     [,42] [,43] [,44] [,45] [,46] [,47] [,48]
ACF  -0.25  0.06 -0.18 -0.03  0.15 -0.22  0.72
PACF  0.06  0.04 -0.07 -0.06  0.02 -0.04  0.10

x = window(hor, start=2002)
par(mfrow = c(2,1))
tsplot(x, main="Hawaiian Quarterly Occupancy Rate", ylab=" % rooms",

ylim=c(62,86), col=gray(.7))
text(x, labels=1:4, col=c(3,4,2,6), cex=.8)

Qx = stl(x,15)$time.series[,1]
tsplot(Qx, main="Seasonal Component", ylab=" % rooms",

ylim=c(-4.7,4.7), col=gray(.7))
text(Qx, labels=1:4, col=c(3,4,2,6), cex=.8)

par(mfrow=c(2,1))
tsplot(cardox, col=4, ylab=expression(CO[2]))
title("Monthly Carbon Dioxide Readings - Mauna Loa Observatory",

cex.main=1)

tsplot(diff(diff(cardox,12)), col=4,

ylab=expression(nabla~nabla[12]~CO[2]))

acf2(diff(diff(cardox,12)))

     [,1]  [,2]  [,3]  [,4]  [,5]  [,6] [,7]  [,8] [,9] [,10] [,11] [,12] [,13] [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21]
ACF  -0.3 -0.03 -0.06  0.00  0.00 -0.01 0.03 -0.01 0.05 -0.03  0.19 -0.47  0.08  0.08  0.00  0.02 -0.05  0.04 -0.04 -0.01 -0.02
PACF -0.3 -0.13 -0.12 -0.07 -0.04 -0.03 0.02  0.01 0.07  0.02  0.22 -0.39 -0.20 -0.03 -0.06 -0.02 -0.07  0.00 -0.01 -0.03  0.02
     [,22] [,23] [,24] [,25] [,26] [,27] [,28] [,29] [,30] [,31] [,32] [,33] [,34] [,35] [,36] [,37] [,38] [,39] [,40] [,41]
ACF   0.02 -0.04  0.00  0.03 -0.06  0.04 -0.03  0.06 -0.03 -0.04  0.04  0.01  0.02  0.00  0.01  0.04 -0.01 -0.01 -0.01 -0.01
PACF  0.00  0.08 -0.25 -0.15 -0.08 -0.04 -0.04 -0.01  0.02 -0.06 -0.02  0.02  0.03  0.07 -0.17 -0.02 -0.06 -0.02 -0.06 -0.03
     [,42] [,43] [,44] [,45] [,46] [,47] [,48]
ACF  -0.03  0.03  0.03 -0.05     0  0.02 -0.02
PACF -0.04 -0.11  0.02 -0.05     0  0.07 -0.14

sarima(cardox, p=0,d=1,q=1, P=0,D=1,Q=1,S=12)

initial  value -0.826338 
iter   2 value -1.059073
iter   3 value -1.093845
iter   4 value -1.116555
iter   5 value -1.124382
iter   6 value -1.126345
iter   7 value -1.127354
iter   8 value -1.127953
iter   9 value -1.127984
iter  10 value -1.127985
iter  10 value -1.127985
iter  10 value -1.127985
final  value -1.127985 
converged
initial  value -1.144615 
iter   2 value -1.148048
iter   3 value -1.148645
iter   4 value -1.149895
iter   5 value -1.150013
iter   6 value -1.150021
iter   7 value -1.150021
iter   7 value -1.150021
iter   7 value -1.150021
final  value -1.150021 
converged
<><><><><><><><><><><><><><>
 
Coefficients: 
     Estimate     SE  t.value p.value
ma1   -0.3869 0.0377 -10.2624       0
sma1  -0.8655 0.0183 -47.2846       0

sigma^2 estimated as 0.0980908 on 766 degrees of freedom 
 
AIC = 0.5456475  AICc = 0.545668  BIC = 0.5637873 
 

sarima(cardox, 1,1,1, 0,1,1,12)

initial  value -0.827261 
iter   2 value -1.034332
iter   3 value -1.066295
iter   4 value -1.094823
iter   5 value -1.108013
iter   6 value -1.115246
iter   7 value -1.116173
iter   8 value -1.120248
iter   9 value -1.120892
iter  10 value -1.121657
iter  11 value -1.122186
iter  12 value -1.124590
iter  13 value -1.125269
iter  14 value -1.125654
iter  15 value -1.125685
iter  16 value -1.125685
iter  17 value -1.125687
iter  18 value -1.125689
iter  19 value -1.125689
iter  20 value -1.125689
iter  20 value -1.125689
iter  20 value -1.125689
final  value -1.125689 
converged
initial  value -1.146682 
iter   2 value -1.150731
iter   3 value -1.152123
iter   4 value -1.152815
iter   5 value -1.153157
iter   6 value -1.153220
iter   7 value -1.153266
iter   8 value -1.153337
iter   9 value -1.153352
iter  10 value -1.153359
iter  11 value -1.153384
iter  12 value -1.153388
iter  13 value -1.153390
iter  14 value -1.153390
iter  14 value -1.153390
iter  14 value -1.153390
final  value -1.153390 
converged
<><><><><><><><><><><><><><>
 
Coefficients: 
     Estimate     SE  t.value p.value
ar1    0.2203 0.0894   2.4660  0.0139
ma1   -0.5797 0.0753  -7.7030  0.0000
sma1  -0.8656 0.0182 -47.5947  0.0000

sigma^2 estimated as 0.09742764 on 765 degrees of freedom 
 
AIC = 0.541514  AICc = 0.5415549  BIC = 0.5657004 
 

sarima.for(cardox, 60, 1,1,1, 0,1,1,12)

$pred
          Jan      Feb      Mar      Apr      May      Jun      Jul      Aug      Sep      Oct      Nov      Dec
2023                            422.5365 423.2516 422.6841 420.8017 418.7844 417.4195 417.6341 419.1994 420.6362
2024 421.8604 422.7144 423.3343 424.7951 425.4936 424.9223 423.0392 421.0216 419.6567 419.8713 421.4367 422.8734
2025 424.0976 424.9517 425.5715 427.0323 427.7308 427.1595 425.2764 423.2588 421.8940 422.1085 423.6739 425.1106
2026 426.3348 427.1889 427.8088 429.2695 429.9680 429.3968 427.5136 425.4961 424.1312 424.3458 425.9111 427.3479
2027 428.5720 429.4261 430.0460 431.5067 432.2052 431.6340 429.7509 427.7333 426.3684 426.5830 428.1483 429.5851
2028 430.8092 431.6633 432.2832                                                                                 

$se
           Jan       Feb       Mar       Apr       May       Jun       Jul       Aug       Sep       Oct       Nov       Dec
2023                               0.3121340 0.3706892 0.4100237 0.4437896 0.4747345 0.5036938 0.5310578 0.5570754 0.5819301
2024 0.6057658 0.6286983 0.6508233 0.6839230 0.7112115 0.7366194 0.7609956 0.7845756 0.8074590 0.8297096 0.8513785 0.8725094
2025 0.8931404 0.9133056 0.9330350 0.9616380 0.9859772 1.0090191 1.0313946 1.0532622 1.0746779 1.0956735 1.1162740 1.1365011
2026 1.1563743 1.1759118 1.1951299 1.2221148 1.2455209 1.2678703 1.2896982 1.3111337 1.3322181 1.3529726 1.3734132 1.3935539
2027 1.4134077 1.4329864 1.4523012 1.4786701 1.5018653 1.5241385 1.5459682 1.5674673 1.5886697 1.6095916 1.6302447 1.6506393
2028 1.6707850 1.6906907 1.7103647                                                                                          

abline(v=2018.9, lty=6)

##-- for comparison, try the first model --##
sarima.for(cardox, 60, 0,1,1, 0,1,1,12) # not shown

$pred
          Jan      Feb      Mar      Apr      May      Jun      Jul      Aug      Sep      Oct      Nov      Dec
2023                            422.4940 423.1899 422.6182 420.7348 418.7172 417.3524 417.5670 419.1324 420.5691
2024 421.7933 422.6474 423.2671 424.7222 425.4182 424.8465 422.9630 420.9455 419.5807 419.7953 421.3606 422.7974
2025 424.0216 424.8757 425.4954 426.9505 427.6464 427.0747 425.1913 423.1738 421.8089 422.0235 423.5889 425.0257
2026 426.2498 427.1039 427.7237 429.1788 429.8747 429.3030 427.4196 425.4020 424.0372 424.2518 425.8172 427.2539
2027 428.4781 429.3322 429.9519 431.4070 432.1030 431.5312 429.6478 427.6303 426.2655 426.4800 428.0454 429.4822
2028 430.7063 431.5605 432.1802                                                                                 

$se
           Jan       Feb       Mar       Apr       May       Jun       Jul       Aug       Sep       Oct       Nov       Dec
2023                               0.3131945 0.3673801 0.4145425 0.4568619 0.4955806 0.5314860 0.5651148 0.5968518 0.6269843
2024 0.6557337 0.6832745 0.7097474 0.7473726 0.7784772 0.8083858 0.8372267 0.8651067 0.8921157 0.9183308 0.9438180 0.9686348
2025 0.9928314 1.0164523 1.0395365 1.0715498 1.0989071 1.1255998 1.1516739 1.1771707 1.2021268 1.2265752 1.2505458 1.2740654
2026 1.2971587 1.3198479 1.3421537 1.3722431 1.3984561 1.4241867 1.4494607 1.4743014 1.4987305 1.5227677 1.5464314 1.5697383
2027 1.5927042 1.6153437 1.6376702 1.6670522 1.6930078 1.7185714 1.7437603 1.7685905 1.7930768 1.8172333 1.8410728 1.8646075
2028 1.8878489 1.9108076 1.9334937                                                                                          

trend = time(cmort); temp = tempr - mean(tempr); temp2 = temp^2
fit = lm(cmort~trend + temp + temp2 + part, na.action=NULL)
acf2(resid(fit), 52) # implies AR2

     [,1] [,2] [,3] [,4]  [,5]  [,6]  [,7]  [,8] [,9] [,10] [,11] [,12] [,13] [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21]
ACF  0.34 0.44 0.28 0.28  0.16  0.12  0.07  0.01 0.03 -0.05 -0.02  0.00 -0.04 -0.02  0.01 -0.05 -0.01 -0.03 -0.06 -0.03 -0.03
PACF 0.34 0.36 0.08 0.06 -0.05 -0.05 -0.02 -0.05 0.02 -0.06  0.00  0.06 -0.02  0.00  0.04 -0.08  0.00 -0.01 -0.06  0.03  0.02
     [,22] [,23] [,24] [,25] [,26] [,27] [,28] [,29] [,30] [,31] [,32] [,33] [,34] [,35] [,36] [,37] [,38] [,39] [,40] [,41]
ACF  -0.05  0.00 -0.02 -0.01 -0.05  0.03 -0.11 -0.02 -0.10 -0.02 -0.09 -0.03 -0.10 -0.08 -0.10 -0.07 -0.07 -0.05 -0.05 -0.04
PACF -0.03  0.05 -0.01  0.00 -0.06  0.06 -0.11  0.00 -0.05  0.05 -0.04  0.04 -0.06 -0.06 -0.05  0.03 -0.02  0.02 -0.01  0.00
     [,42] [,43] [,44] [,45] [,46] [,47] [,48] [,49] [,50] [,51] [,52]
ACF  -0.03 -0.04  0.05  0.00  0.04  0.08  0.07  0.06  0.05  0.07  0.02
PACF  0.00 -0.01  0.08 -0.01 -0.01  0.08 -0.01 -0.01 -0.03  0.03 -0.04

sarima(cmort, 2,0,0, xreg=cbind(trend, temp, temp2, part) )

initial  value 1.849900 
iter   2 value 1.733730
iter   3 value 1.701063
iter   4 value 1.682463
iter   5 value 1.657377
iter   6 value 1.652444
iter   7 value 1.641726
iter   8 value 1.635302
iter   9 value 1.630848
iter  10 value 1.629286
iter  11 value 1.628731
iter  12 value 1.628646
iter  13 value 1.628634
iter  14 value 1.628633
iter  15 value 1.628632
iter  16 value 1.628628
iter  17 value 1.628627
iter  18 value 1.628627
iter  19 value 1.628626
iter  20 value 1.628625
iter  21 value 1.628622
iter  22 value 1.628618
iter  23 value 1.628614
iter  24 value 1.628612
iter  25 value 1.628611
iter  26 value 1.628610
iter  27 value 1.628610
iter  28 value 1.628609
iter  29 value 1.628609
iter  30 value 1.628608
iter  31 value 1.628608
iter  32 value 1.628608
iter  32 value 1.628608
iter  32 value 1.628608
final  value 1.628608 
converged
initial  value 1.630401 
iter   2 value 1.630393
iter   3 value 1.630382
iter   4 value 1.630381
iter   5 value 1.630375
iter   6 value 1.630370
iter   7 value 1.630362
iter   8 value 1.630354
iter   9 value 1.630349
iter  10 value 1.630347
iter  11 value 1.630347
iter  12 value 1.630347
iter  13 value 1.630347
iter  14 value 1.630346
iter  15 value 1.630346
iter  16 value 1.630346
iter  17 value 1.630346
iter  17 value 1.630346
iter  17 value 1.630346
final  value 1.630346 
converged
<><><><><><><><><><><><><><>
 
Coefficients: 
           Estimate       SE t.value p.value
ar1          0.3848   0.0436  8.8329  0.0000
ar2          0.4326   0.0400 10.8062  0.0000
intercept 3075.1482 834.7278  3.6840  0.0003
trend       -1.5165   0.4226 -3.5881  0.0004
temp        -0.0190   0.0495 -0.3837  0.7014
temp2        0.0154   0.0020  7.6117  0.0000
part         0.1545   0.0272  5.6803  0.0000

sigma^2 estimated as 26.01476 on 501 degrees of freedom 
 
AIC = 6.130066  AICc = 6.130507  BIC = 6.196687 
 

library(zoo)
lag2.plot(Hare, Lynx, 5) # lead-lag relationship

pp = as.zoo(ts.intersect(Lynx, HareL1 = lag(Hare,-1)))

summary(reg <- lm(pp$Lynx~ pp$HareL1)) # results not displayed

Call:
lm(formula = pp$Lynx ~ pp$HareL1)

Residuals:
    Min      1Q  Median      3Q     Max 
-31.107 -11.511  -2.193   7.579  45.632 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 15.84712    2.75817   5.746 1.29e-07 ***
pp$HareL1    0.27265    0.04727   5.768 1.17e-07 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 16.25 on 88 degrees of freedom
Multiple R-squared:  0.2744,    Adjusted R-squared:  0.2661 
F-statistic: 33.27 on 1 and 88 DF,  p-value: 1.172e-07

acf2(resid(reg)) # in Figure 5.19

     [,1]  [,2]  [,3]  [,4]  [,5]  [,6]  [,7]  [,8] [,9] [,10] [,11] [,12] [,13] [,14] [,15] [,16] [,17] [,18] [,19] [,20]
ACF  0.73  0.30 -0.10 -0.35 -0.48 -0.43 -0.21  0.03 0.21  0.28  0.27  0.09 -0.16 -0.32 -0.35 -0.26 -0.10  0.14  0.34  0.37
PACF 0.73 -0.52 -0.18 -0.05 -0.25  0.08  0.13 -0.08 0.04  0.03  0.00 -0.24 -0.10  0.07 -0.13  0.04  0.06  0.08  0.05 -0.10

( reg2 = sarima(pp$Lynx, 2,0,0, xreg=pp$HareL1 ))

initial  value 2.758050 
iter   2 value 2.608114
iter   3 value 2.338840
iter   4 value 2.248812
iter   5 value 2.046836
iter   6 value 2.045593
iter   7 value 2.045218
iter   8 value 2.045114
iter   9 value 2.045060
iter  10 value 2.045055
iter  11 value 2.045055
iter  12 value 2.045055
iter  13 value 2.045055
iter  13 value 2.045055
iter  13 value 2.045055
final  value 2.045055 
converged
initial  value 2.056766 
iter   2 value 2.056665
iter   3 value 2.056583
iter   4 value 2.056581
iter   5 value 2.056580
iter   6 value 2.056580
iter   6 value 2.056580
iter   6 value 2.056580
final  value 2.056580 
converged
<><><><><><><><><><><><><><>
 
Coefficients: 
          Estimate     SE t.value p.value
ar1         1.3258 0.0732 18.1184  0.0000
ar2        -0.7143 0.0731 -9.7689  0.0000
intercept  25.1319 2.5469  9.8676  0.0000
xreg        0.0692 0.0318  2.1727  0.0326

sigma^2 estimated as 59.57091 on 86 degrees of freedom 
 
AIC = 7.062148  AICc = 7.067377  BIC = 7.201026 
 
$fit

Call:
arima(x = xdata, order = c(p, d, q), seasonal = list(order = c(P, D, Q), period = S), 
    xreg = xreg, transform.pars = trans, fixed = fixed, optim.control = list(trace = trc, 
        REPORT = 1, reltol = tol))

Coefficients:
         ar1      ar2  intercept    xreg
      1.3258  -0.7143    25.1319  0.0692
s.e.  0.0732   0.0731     2.5469  0.0318

sigma^2 estimated as 59.57:  log likelihood = -312.8,  aic = 635.59

$degrees_of_freedom
[1] 86

$ttable
          Estimate     SE t.value p.value
ar1         1.3258 0.0732 18.1184  0.0000
ar2        -0.7143 0.0731 -9.7689  0.0000
intercept  25.1319 2.5469  9.8676  0.0000
xreg        0.0692 0.0318  2.1727  0.0326

$ICs
     AIC     AICc      BIC 
7.062148 7.067377 7.201026 

prd = Lynx - resid(reg2$fit) # prediction (resid = obs - pred)
prde = sqrt(reg2$fit$sigma2) # prediction error
tsplot(prd, lwd=2, col=rgb(0,0,.9,.5), ylim=c(-20,90), ylab="Lynx")
points(Lynx, pch=16, col=rgb(.8,.3,0))

x = time(Lynx)[-1]
xx = c(x, rev(x))
yy = c(prd - 2*prde, rev(prd + 2*prde))
polygon(xx, yy, border=8, col=rgb(.4, .5, .6, .15))
mtext(expression(""%*% 10^3), side=2, line=1.5, adj=.975)

legend("topright", legend=c("Predicted", "Observed"), lty=c(1,NA),
lwd=2, pch=c(NA,16), col=c(4,rgb(.8,.3,0)), cex=.9)

t = seq(0, 24, by=.01)
X = cos(2*pi*t*1/2) # 1 cycle every 2 hours
tsplot(t, X, xlab="Hours")
T = seq(1, length(t), by=250) # observed every 2.5 hrs
points(t[T], X[T], pch=19, col=4)

lines(t, cos(2*pi*t/10), col=4)
axis(1, at=t[T], labels=FALSE, lwd.ticks=2, col.ticks=2)

abline(v=t[T], col=rgb(1,0,0,.2), lty=2)

x1 = 2*cos(2*pi*1:100*6/100) + 3*sin(2*pi*1:100*6/100)
x2 = 4*cos(2*pi*1:100*10/100) + 5*sin(2*pi*1:100*10/100)
x3 = 6*cos(2*pi*1:100*40/100) + 7*sin(2*pi*1:100*40/100)
x = x1 + x2 + x3
par(mfrow=c(2,2))
tsplot(x1, ylim=c(-10,10), main=expression(omega==6/100~~~A^2==13))
tsplot(x2, ylim=c(-10,10), main=expression(omega==10/100~~~A^2==41))

tsplot(x3, ylim=c(-10,10), main=expression(omega==40/100~~~A^2==85))
tsplot(x, ylim=c(-16,16), main="sum")

P = Mod(fft(x)/sqrt(100))^2 # periodogram
sP = (4/100)*P # scaled peridogram
Fr = 0:99/100 # fundamental frequencies
tsplot(Fr, sP, type="o", xlab="frequency", ylab="scaled periodogram",

col=4, ylim=c(0,90))

abline(v=.5, lty=5)

abline(v=c(.1,.3,.7,.9), lty=1, col=gray(.9))
axis(side=1, at=seq(.1,.9,by=.2))

par(mfrow=c(3,2), mar=c(1.5,2,1,0)+1, mgp=c(1.6,.6,0))
for(i in 4:9){
mvspec(fmri1[,i], main=colnames(fmri1)[i], ylim=c(0,3), xlim=c(0,.2),

col=rgb(.05,.6,.75), lwd=2, type="o", pch=20)
abline(v=1/32, col="dodgerblue", lty=5) # stimulus frequency
}

par(mfrow=c(3,1))
arma.spec(main="White Noise", col=4)
arma.spec(ma=.5, main="Moving Average", col=4)

arma.spec(ar=c(1,-.9), main="Autoregression", col=4)

par(mfrow=c(3,1))
tsplot(soi, col=4, main="SOI")
tsplot(diff(soi), col=4, main="First Difference")

k = kernel("modified.daniell", 6) # MA weights
tsplot(kernapply(soi, k), col=4, main="Seasonal Moving Average")
##-- frequency responses --##
par(mfrow=c(2,1))

w = seq(0, .5, by=.01)
FRdiff = abs(1-exp(2i*pi*w))^2
tsplot(w, FRdiff, xlab="frequency", main="High Pass Filter")
u = cos(2*pi*w)+cos(4*pi*w)+cos(6*pi*w)+cos(8*pi*w)+cos(10*pi*w)
FRma = ((1 + cos(12*pi*w) + 2*u)/12)^2
tsplot(w, FRma, xlab="frequency", main="Low Pass Filter")

library(fGarch)
library(tseries)
library(arfima)
library(tsDyn)
library(ggplot2)
library(dplyr)

######################### Warning from 'xts' package ##########################
#                                                                             #
# The dplyr lag() function breaks how base R's lag() function is supposed to  #
# work, which breaks lag(my_xts). Calls to lag(my_xts) that you type or       #
# source() into this session won't work correctly.                            #
#                                                                             #
# Use stats::lag() to make sure you're not using dplyr::lag(), or you can add #
# conflictRules('dplyr', exclude = 'lag') to your .Rprofile to stop           #
# dplyr from breaking base R's lag() function.                                #
#                                                                             #
# Code in packages is not affected. It's protected by R's namespace mechanism #
# Set `options(xts.warn_dplyr_breaks_lag = FALSE)` to suppress this warning.  #
#                                                                             #
###############################################################################

Attaching package: ‘dplyr’

The following objects are masked from ‘package:xts’:

    first, last

The following objects are masked from ‘package:stats’:

    filter, lag

The following objects are masked from ‘package:base’:

    intersect, setdiff, setequal, union

renv::init()

renv: Project Environments for R

Welcome to renv! It looks like this is your first time using renv.
This is a one-time message, briefly describing some of renv's functionality.

renv will write to files within the active project folder, including:

  - A folder 'renv' in the project directory, and
  - A lockfile called 'renv.lock' in the project directory.

In particular, projects using renv will normally use a private, per-project
R library, in which new packages will be installed. This project library is
isolated from other R libraries on your system.

In addition, renv will update files within your project directory, including:

  - .gitignore
  - .Rbuildignore
  - .Rprofile

Finally, renv maintains a local cache of data on the filesystem, located at:

  - "~/.cache/R/renv"

This path can be customized: please see the documentation in `?renv::paths`.

Please read the introduction vignette with `vignette("renv")` for more information.
You can browse the package documentation online at https://rstudio.github.io/renv/.

y

- "~/.cache/R/renv" has been created.
- Linking packages into the project library ... [73/128] [74/128] [75/128] [76/128] [77/128] [78/128] [79/128] [80/128] [81/128] [82/128] [83/128] [84/128] [85/128] [86/128] [87/128] [88/128] [89/128] [90/128] [91/128] [92/128] [93/128] [94/128] [95/128] [96/128] [97/128] [98/128] [99/128] [100/128] [101/128] [102/128] [103/128] [104/128] [105/128] [106/128] [107/128] [108/128] [109/128] [110/128] [111/128] [112/128] [113/128] [114/128] [115/128] [116/128] [117/128] [118/128] [119/128] [120/128] [121/128] [122/128] [123/128] [124/128] [125/128] [126/128] [127/128] [128/128] Done!
- Resolving missing dependencies ... 
# Installing packages --------------------------------------------------------
The following package(s) will be updated in the lockfile:

# CRAN -----------------------------------------------------------------------
- abind           [* -> 1.4-5]
- arfima          [* -> 1.8-1]
- astsa           [* -> 2.1]
- backports       [* -> 1.4.1]
- base64enc       [* -> 0.1-3]
- boot            [* -> 1.3-28]
- brio            [* -> 1.1.5]
- broom           [* -> 1.0.5]
- bslib           [* -> 0.7.0]
- cachem          [* -> 1.0.8]
- callr           [* -> 3.7.6]
- car             [* -> 3.1-2]
- carData         [* -> 3.0-5]
- cli             [* -> 3.6.2]
- codetools       [* -> 0.2-18]
- colorspace      [* -> 2.1-0]
- cpp11           [* -> 0.4.7]
- crayon          [* -> 1.5.2]
- curl            [* -> 5.2.1]
- cvar            [* -> 0.5]
- desc            [* -> 1.4.3]
- deSolve         [* -> 1.40]
- diffobj         [* -> 0.3.5]
- digest          [* -> 0.6.35]
- dplyr           [* -> 1.1.4]
- dynlm           [* -> 0.3-6]
- ellipsis        [* -> 0.3.2]
- evaluate        [* -> 0.23]
- fansi           [* -> 1.0.6]
- farver          [* -> 2.1.1]
- fastICA         [* -> 1.2-4]
- fastmap         [* -> 1.1.1]
- fBasics         [* -> 4032.96]
- fGarch          [* -> 4033.92]
- fontawesome     [* -> 0.5.2]
- foreach         [* -> 1.5.2]
- forecast        [* -> 8.22.0]
- fracdiff        [* -> 1.5-3]
- fs              [* -> 1.6.4]
- gbutils         [* -> 0.5]
- generics        [* -> 0.1.3]
- ggplot2         [* -> 3.5.1]
- glue            [* -> 1.7.0]
- gss             [* -> 2.2-7]
- gtable          [* -> 0.3.5]
- highr           [* -> 0.10]
- htmltools       [* -> 0.5.8.1]
- isoband         [* -> 0.2.7]
- iterators       [* -> 1.0.14]
- jquerylib       [* -> 0.1.4]
- jsonlite        [* -> 1.8.8]
- knitr           [* -> 1.46]
- labeling        [* -> 0.4.3]
- lattice         [* -> 0.20-45]
- lifecycle       [* -> 1.0.4]
- lme4            [* -> 1.1-35.3]
- lmtest          [* -> 0.9-40]
- ltsa            [* -> 1.4.6]
- magrittr        [* -> 2.0.3]
- MASS            [* -> 7.3-55]
- Matrix          [* -> 1.6-5]
- MatrixModels    [* -> 0.5-3]
- memoise         [* -> 2.0.1]
- mgcv            [* -> 1.8-39]
- mime            [* -> 0.12]
- minqa           [* -> 1.2.6]
- mnormt          [* -> 2.1.1]
- munsell         [* -> 0.5.1]
- nlme            [* -> 3.1-155]
- nloptr          [* -> 2.0.3]
- nnet            [* -> 7.3-17]
- numDeriv        [* -> 2016.8-1.1]
- pbkrtest        [* -> 0.5.2]
- pillar          [* -> 1.9.0]
- pkgbuild        [* -> 1.4.4]
- pkgconfig       [* -> 2.0.3]
- pkgload         [* -> 1.3.4]
- praise          [* -> 1.0.0]
- processx        [* -> 3.8.4]
- ps              [* -> 1.7.6]
- purrr           [* -> 1.0.2]
- quadprog        [* -> 1.5-8]
- quantmod        [* -> 0.4.26]
- quantreg        [* -> 5.97]
- R6              [* -> 2.5.1]
- rappdirs        [* -> 0.3.3]
- rbibutils       [* -> 2.2.16]
- RColorBrewer    [* -> 1.1-3]
- Rcpp            [* -> 1.0.12]
- RcppArmadillo   [* -> 0.12.8.2.1]
- RcppEigen       [* -> 0.3.4.0.0]
- Rdpack          [* -> 2.6]
- rematch2        [* -> 2.1.2]
- renv            [* -> 1.0.7]
- rlang           [* -> 1.1.4]
- rmarkdown       [* -> 2.26]
- rprojroot       [* -> 2.0.4]
- sandwich        [* -> 3.1-0]
- sass            [* -> 0.4.9]
- scales          [* -> 1.3.0]
- SparseM         [* -> 1.81]
- spatial         [* -> 7.3-15]
- stabledist      [* -> 0.7-1]
- stringi         [* -> 1.8.3]
- stringr         [* -> 1.5.1]
- strucchange     [* -> 1.5-3]
- survival        [* -> 3.2-13]
- testthat        [* -> 3.2.1.1]
- tibble          [* -> 3.2.1]
- tidyr           [* -> 1.3.1]
- tidyselect      [* -> 1.2.1]
- timeDate        [* -> 4032.109]
- timeSeries      [* -> 4032.109]
- tinytex         [* -> 0.50]
- tsDyn           [* -> 11.0.4.1]
- tseries         [* -> 0.10-56]
- tseriesChaos    [* -> 0.1-13.1]
- TTR             [* -> 0.24.4]
- urca            [* -> 1.3-4]
- utf8            [* -> 1.2.4]
- vars            [* -> 1.6-1]
- vctrs           [* -> 0.6.5]
- viridisLite     [* -> 0.4.2]
- waldo           [* -> 0.5.2]
- withr           [* -> 3.0.0]
- xfun            [* -> 0.43]
- xts             [* -> 0.14.0]
- yaml            [* -> 2.3.8]
- zoo             [* -> 1.8-12]

The version of R recorded in the lockfile will be updated:
- R               [* -> 4.1.2]

- Lockfile written to "~/water-quality-time-series/renv.lock".